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Implicitly  Coupling  Heat  Conduction  into  a 

Zone  Fire  Model 


William  F.  Moss*  Glenn  P.  Forney* 


Abstract 

This  report  examines  several  methods  for  coupling  the  partial  differential 
equations  that  arise  in  conductive  heat  transfer  with  the  ordinary  differential 
equations  that  arise  in  zone  fire  modeling.  Two  existing  algorithms  (method 
of  lines  and  time  splitting)  are  discussed  and  a new  strategy  is  proposed  for 
performing  this  coupling.  This  strategy  couples  the  wall  surface  temperature 
rather  than  the  entire  wall  temperature  profile  with  the  other  zone  fire  modeling 
solution  variables  by  requiring  that  the  wall  surface  temperature  gradient  and 
the  incident  heat  flux  (sum  of  convective  and  net  radiative  flux)  satisfy  Fourier’s 
law,  q"  = —Kdu/dx. 

Two  prototype  Are  models  were  written  to  test  the  ideas  discussed  in  this 
report.  The  first,  CONRADl,  implements  the  method  of  lines  strategy  for 
solving  heat  conduction.  The  second,  C0NRAD2,  implements  the  new  strategy. 
Though  inefficient,  CONRADl  uses  well  established  numerical  techniques  and 
therefore  serves  as  a benchmark  to  test  the  numerical  ideas  implemented  in 
C0NRAD2.  Both  programs  use  the  stiff  differential-algebraic  equation  solver 
DASSL.  Supporting  numerical  results  are  presented. 


1 Introduction 

In  a zone  fire  model,  each  room  in  a building  is  divided  into  a relatively  smoke-laden 
upper  layer  and  a relatively  clear  lower  layer.  Temperatures  in  the  ceiling,  floor,  and 
vertical  walls  of  a room  must  be  computed  in  order  to  adequately  account  for  heat 
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exchange  due  to  convective  and  radiative  heat  transfer.  The  ceiling,  floor,  and  vertical 
walls  of  each  room  are  usually  divided  into  “wall”  segments.  Here  the  term  “wall” 
may  mean  the  ceiling,  the  floor,  or  the  vertical  walls  of  a room.  Wall  temperatures  are 
usually  computed  by  solving  a series  of  1-d  heat  conduction  problems;  one  problem 
per  wall  segment.  The  gas  layers  in  a zone  fire  model  are  often  modeled  using  ordinary 
dilferential  equations  (ODE)  which  are  derived  from  conservation  of  mass  and  energy. 
CCFM. VENTS  [1]  and  FAST  (renamed  CFAST)  [2]  are  two  examples  of  zone  fire 
models.  CCFM. VENTS  uses  pressure,  layer  interface  height,  and  upper  and  lower 
layer  masses  as  solution  variables,  while  CFAST  uses  pressure,  upper  layer  volume, 
and  upper  and  lower  layer  temperatures.  Reference  [3]  shows  how  these  two  and 
many  other  formulations  are  equivalent  in  the  sense  that  one  can  be  converted  into 
another  using  the  ideal  gas  law  and  definitions  of  physical  properties  such  as  internal 
energy  and  density. 

The  heat  conduction  problems  for  the  wall  segments  are  formulated  using  the  heat 
conduction  equation,  a partial  differential  equation  (PDF)  [4],  which  can  be  derived 
using  the  conservation  of  energy.  The  gas  solution  variables  and  the  wall  segment 
temperature  profiles  are  presumed  known  at  some  time  t.  The  ODE’s  and  PDF’s  are 
used  to  advance  the  solution  variables  to  some  later  time  i + Ai.  The  gas  layers  and 
the  wall  segments  or  equivalently  the  ODE’s  and  the  PDF’s  are  coupled  via  convective 
and  radiative  heat  transfer  terms. 

Procedures  for  solving  1-d  heat  conduction  problems  are  well  known.  For  finite 
difference  methods  such  as  backward  difference  (fully  implicit),  forward  difference 
(fully  explicit)  or  Crank-Nicolson  see  [5].  For  finite  element  methods  see  [6]. 

The  question  addressed  by  this  report  is  how  to  couple  the  ODE’s  from  the  gas 
layers  with  the  PDF’s  from  the  wall  segments  in  a numerically  accurate,  robust  and 
efficient  manner. 

Two  prototype  fire  codes  were  written  to  test  the  ideas  presented  in  this  report. 
CONRAD  1 implements  the  method  of  lines  strategy  for  solving  the  heat  equation 
using  standard  cubic  Hermite  polynomials  to  represent  the  unknown  wall  temper- 
ature profile.  C0NRAD2  implements  a new  method  called  “gradient  matching” 
which  is  based  on  an  implicitly  defined  functional  equation  approach.  Though  in- 
efficient, CONRAD  1 uses  well  established  numerical  techniques  and  therefore  serves 
as  a benchmark  to  test  the  new  numerical  ideas  implemented  in  C0NRAD2.  CON- 
RAD 1 and  CONRAD2  both  use  the  stiff  differential-algebraic  equation  (DAE)  solver 
DASSL[7,  8].  These  zone  fire  models,  documented  in  [9],  use  the  same  solution  vari- 
ables as  CCFM. VENTS  which  are  denoted  by  P,  y,  mu  and  mu.  Here  these  variables 
are  referred  to  as  the  gas  solution  variables.  The  procedures  discussed  in  this  report 
for  coupling  the  gas  zone  properties  with  the  wall  temperatures  will  work  for  many 
other  formulations  of  the  gas  solution  variables. 
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2 Background 

Two  standard  approaches  for  coupling  zone  fire  modeling  differential  equations  with 
1-d  heat  conduction  problems  are  time  splitting  and  the  method  of  lines  (MOL). 
These  two  methods  will  be  briefly  discussed  in  the  next  two  subsections. 


2.1  Time  Splitting 

Time  splitting  makes  the  assumption  that  two  or  more  phenomena  change  over  sig- 
nificantly different  time  scales.  For  example,  in  the  zone  fire  modeling  case,  it  can 
often  be  assumed  that  the  characteristic  time  scale  for  wall  segment  temperature  pro- 
files is  much  longer  than  that  for  the  gas  solution  variables.  Suppose  that  the  gas 
solution  variables  and  the  wall  segment  temperature  profiles  are  known  at  time  t.  If 
the  characteristic  time  scale  for  wall  segment  temperature  profiles  is  At,  then  wall 
segment  temperature  profiles  would  be  solved  over  the  time  interval  {t,t  + At).  This 
time  interval  would  then  be  further  subdivided  in  order  to  solve  for  the  shorter  time 
scale  phenomena.  The  longer  time  interval  is  often  called  the  outer  time  step  and  the 
shorter  interval  is  called  the  inner  time  step. 

Referring  to  Figure  1,  time  splitting  advances  the  short  time  scale  phenomena 
from  time  t to  time  t-1- At  in  a series  of  time  steps  chosen  sufficiently  small  by  the  solver 
to  satisfy  the  error  criteria.  Next  the  wall  segment  temperature  profiles,  are  advanced 
from  time  t to  time  t -f-  At.  The  outer  time  stepsize  At  must  be  chosen  sufficiently 
small  so  that  the  computed  wall  segment  temperature  profiles  are  consistent  with  the 
gas  solution  variables  at  time  t + At.  By  consistent  it  is  meant  that  the  heat  flux 
striking  each  wall  segment,  q",  is  related  to  the  temperature  gradient,  dujdx,  at  the 
surface  of  the  wall  segment  via  Fourier’s  law,  q"  = —Kdufdx,  where  K is  the  thermal 
conductivity  of  the  wall. 

The  wall  segment  temperature  profiles  are  advanced  using  a flux  boundary  con- 
dition. The  flux  at  the  interior  surface  of  a wall  segment  is  the  sum  of  convective 
and  net  radiative  heat  fluxes.  The  question  then  is  what  flux  to  use  to  advance  the 
solution;  the  flux  at  time  t,  the  flux  at  time  t + At,  or  some  combination.  The  flux 
at  time  t At  is  not  known  until  both  the  inner  variables  {P,  y,  mu  and  m^)  and 
the  outer  variables  (wall  segment  temperature  profiles)  are  determined  ai  t At. 
Consequently,  an  iterative  procedure  must  be  used  to  insure  that  the  flux  striking  the 
interior  surface  of  a wall  segment  is  consistent  with  the  temperature  gradient  there. 

The  method  of  time  splitting  does  not  work  well  when  the  time  scales  are  close 
which  can  occur  when  wall  materials  are  thin  and/or  highly  conductive.  Time  split- 
ting is  also  difficult  to  implement  efficiently  since  it  is  not  clear  what  time  stepsizes 
should  be  used.  A time  stepsize  chosen  too  small  will  result  in  inefficiency  and  time 
stepsize  chosen  too  large  will  result  in  unnecessarily  inaccurate  answers.  The  former 
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Inner  Solution  Variables 
Pressure,  layer  interface 
lower/upper  layer  mass 


Figure  1:  A General  Flowchart  for  the  Method  of  Time  Splitting 


can  easily  occur  as  a fire  simulation  approaches  steady  state.  In  this  case,  solution 
variables  do  not  change  much  and  inefficiency  occurs  because  of  restrictive  time  step- 
size  selections. 

2.2  Method  of  Lines 

The  MOL  consists  of  converting  the  heat  conduction  partial  differential  equation  into 
a system  of  ordinary  differential  equations.  The  unknown  wall  segment  temperature 
profile  is  expanded  as 

N 

u{x,t)  = J2^k{t)Bk{x)  (1) 

k=l 

where  u(x,  t)  is  an  approximation  to  the  unknown  temperature  profile  at  a distance  x 
into  the  wall  segment  at  time  t.  The  functions  Bk{x)  are  known  basis  functions  and 
the  ak{t)  are  unknown  coefficient  functions.  Equation  (1)  is  substituted  into  the  heat 
equation 

du{x,t)  d^u{x,t) 
dt  ^ 5x2  ’ 

where  a is  the  thermal  diffusivity,  to  obtain  a differential  equation  for  the  unknown 
coefficients  a^: 

^ a'f^{t)Bk{x)  = aY^  ak{t)B',^{x)  . 

A:=l  k=l 

A system  of  ODE’s  is  obtained  by  requiring  that  the  above  equation  be  satisfied 
at  a set  of  collocation  points  x = Xi, . . . , XAr-2-  The  flux  boundary  conditions  at 
each  end  of  the  wall  segment  generate  two  additional  equations.  Many  variations  of 
the  MOL  algorithm  can  be  derived  by  choosing  different  basis  functions  (B-Splines, 
cubic  Hermite  interpolating  polynomials,  trigonometric  polynomials,  for  example) 
and  different  sets  of  collocation  points.  The  MOL  algorithm  in  CONRAD  1 uses 
the  standard  basis  functions  for  cubic  Hermite  polynomial  interpolation,  and  it  uses 
Gaussian  points  for  collocation  points.  This  is  discussed  in  Appendix  A. 

In  CONRADl,  the  ODE’s  for  the  the  gas  solution  variables  and  the  coefficients 
ak  are  solved  simultaneously.  Heat  conduction  and  gas  phenomena  are  coupled  via 
radiative  and  convective  heat  transfer. 

In  CONRADl  and  C0NRAD2  the  four  gas  solution  variables  have  four  associated 
ODE’s.  In  CONRADl,  the  number  of  additional  ODE’s  required  per  room  is  NM 
where  N is  the  number  of  basis  functions  used  to  represent  the  wall  segment  temper- 
ature profiles  and  M is  the  number  of  wall  segments  per  room.  Even  for  moderate 


5 


values  of  M,  TV,  and  the  number  of  rooms,  the  run-times  required  for  CONRADl 
can  easily  increase  by  an  order  of  magnitude  over  a simulation  with  the  conduction 
submodel  turned  off. 

CONRADl  advances  the  gas  solution  variables  along  with  the  wall  temperature 
profiles  from  time  t to  t At.  The  time  stepsize.  At,  is  chosen  by  the  DAE  solver 
to  satisfy  appropriate  error  criteria.  The  flux  boundary  conditions  form  part  of  the 
system  of  equations  being  solved.  As  a result,  the  inconsistencies  that  could  arise 
using  time  splitting  will  not  occur.  This  gain  in  consistency,  however,  is  at  the 
expense  of  increased  computational  requirements. 

As  is  the  case  here,  the  MOL  often  requires  a DAE  solver  because  the  boundary 
conditions  are  algebraic  equations  which  cannot  be  easily  converted  into  ODE’s.  For 
simple  problems  in  which  the  boundary  conditions  can  be  easily  differentiated  with 
respect  to  time,  a standard  (stiff)  ODE  solver  can  be  used. 

3 Coupling  Heat  Conduction  Using  Gradient  Match- 
ing 

This  report  presents  a new  strategy  for  coupling  1-d  heat  conduction  problem  with 
the  ODE’s  for  the  gas  solution  variables.  Only  the  temperature  of  the  interior  wall 
segment  surface  (as  opposed  to  the  entire  wall  segment  temperature  profile)  is  di- 
rectly coupled  with  the  gas  layers  and  fire  through  convective  and  radiative  heat 
transfer.  This  observation  has  been  exploited  to  design  a new  and  efficient  algorithm 
for  coupling  heat  conduction  with  a zone  fire  model. 

3.1  Description 

The  method  for  modeling  heat  conduction  discussed  in  this  section  couples  the  wall 
segment  surface  temperatures,  rather  than  the  entire  wall  segment  temperature  pro- 
file, with  the  gas  solution  variables  by  requiring  that  the  wall  segment  surface  temper- 
ature gradient,  and  the  incident  heat  flux  (sum  of  convective  and  net  radiative 

flux),  q”  satisfy  Fourier’s  law 

//  _ T^du{x,t) 

^ dx 

at  the  wall  boundaries  x = 0 and  x = W where  K is  the  thermal  conductivity  of 
the  wall  material  and  W is  the  wall  thickness.  This  solution  strategy  requires  a 
DAE  solver  that  can  simultaneously  solve  both  differential  (gas  ODE’s)  and  algebraic 
equations  (Fourier’s  law).  With  this  method,  only  one  or  two  extra  equations  are 
required  per  wall  segment  (two  if  both  the  interior  and  exterior  wall  segment  surface 
temperatures  are  computed).  This  solution  strategy  is  computationally  more  efficient 
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Figure  2:  A General  Flowchart  for  the  Method  of  Gradient  Matching 


than  the  method  of  lines  since  fewer  equations  need  to  be  solved.  Wall  segment 
temperature  profiles,  however,  still  have  to  be  stored  so  there  is  no  decrease  in  storage 
requirements. 

Consider  a room  with  a single  wall  segment  with  both  interior  and  exterior  wall 
segment  surface  temperatures  computed.  In  this  case,  there  will  be  six  solution 
variables,  the  four  gas  solution  variables,  P,  y,  mt/,  m^,  plus  the  two  wall  segment 
surface  temperatures.  There  will  be  six  equations  to  solve,  four  ODE’s  associated 
with  the  gas  solution  variables,  plus  two  algebraic  equations  consisting  of  Fourier’s 
law  applied  at  the  surfaces  of  the  wall  segment.  Referring  to  Figure  2,  the  gradient 
matching  method  assumes  that  the  gas  solution  variables  and  the  wall  segment 
temperature  profile  are  known  at  time  t.  The  DAE  solver  will  make  an  initial  guess 
at  values  for  the  solution  variables  at  time  t + At.  Based  on  the  wall  segment  surface 
temperatures  at  times  t and  < -f  At,  the  wall  segment  temperature  profile  is  advanced 
to  time  t + At.  Next,  the  wall  segment  surface  temperature  gradients  are  estimated  at 
time  t + At.  Finally,  the  residuals  are  computed  at  time  t + At  for  the  six  equations 
including  the  two  Fourier  law  equations.  The  DAE  solver  adjusts  the  stepsize  At 
until  the  residuals  for  all  six  equations  are  below  an  error  tolerance. 
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3.2  Theoretical  Justification  for  the  Algorithm 

The  basic  idea  is  to  transform  the  initial-boundary  value  problem  for  each  wall  seg- 
ment into  a pair  of  functional  equations.  To  explain  further,  consider  the  following 
standard  problem.  Find  the  temperature  profile  u{x,t)  that  satisfies 


du{x^  t) 
dt 

d^uix^t) 

= a — — — , 0<x<IF,f>0 
ox^ 

(3) 

u(x,  0) 

= g{x)^  0 < X < VF 

(4) 

u(0, t) 

= fo{t),  t > 0 

(5) 

u{W,t) 

= i > 0, 

(6) 

where  g,  /o,  and  /i  are  given  continuous  functions,  a = — is  the  thermal  diffusivity, 
K is  the  thermal  conductivity,  p is  the  density  and  c is  the  specific  heat  of  the  wall 
material.  The  existence  and  uniqueness  theory  (see  [10])  for  this  problem  shows  that 
u{x^t)  is  uniquely  determined  by  the  initial  temperature  profile,  ^r,  and  the  interior 
and  exterior  temperature  boundary  functions  /oC"?")  and  /i(t)  for  0 < r < f.  In  other 
words,  there  is  a functional  relation  between  the  initial  temperature  profile,  boundary 
functions,  and  the  solution  of  the  form 

u(x,<)  = G'[^,/o(0  < T < f),/i(0  < T < f)](x,f)  . (7) 

In  [10]  the  Laplace  transform  is  used  to  construct  the  functional  G.  In  the  special 
case  when  the  initial  temperature  is  constant,  say  u(x,0)  = Tamb,  the  functional  G 
can  be  written  in  terms  of  convolution  integrals.  With  a = VF  = 1,  the  solution  has 
the  form 

U{x,t)  = Tamb  - ^ /o(f  - ~ ’ 

where  03(a;,f)  is  a theta  function  [11]  which  can  be  expressed  as 

OO 

9^{x^  f)  = 1 -f  2 ^ cos  2Tvkx  . 

k=l 


From  equation  (7),  substituting  a:  = 0 for  the  interior  and  x = W for  the  exterior 
wall  segment  boundaries,  it  follows  that 


1^(0, f)  = ^G[gjo{0  <T  <T  <t)]{0,t) 

= ^G'[s,/o(0<T<i),/.(0<T<i)](Vl/,O. 
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For  each  wall  segment  in  a room,  two  functional  equations  are  added  to  four  gas 
ode’s.  These  equations,  based  on  Fourier’s  law  (equation  (2)),  have  the  form 

A"— GfTamb,  w(0,0  < r < <),u(VF,0  < r < t)]{0,t) 

+ 9"(u(0,  t),  mL{t),  mu{t),  y{t),  P{t))  = 0 (8) 

K—G[T^^,u{0,0  <T<  t),u{W,0  <T<  t)]{W,t) 

+ q'\u{W,  t),  mL{t),  mu{t),  y{t),  P{t))  = 0 . (9) 

The  expressions  u(0,0  < t < t)  and  u(VF,0  < t < f)  in  equations  (8)  and  (9)  denote 
the  temperatures  of  the  interior  and  exterior  wall  segment  surfaces  over  the  time 
interval  (0,i)  which  makes  these  functional  equations. 

Returning  to  the  one  room,  one  wall  segment  discussion  of  the  previous  section, 
there  are  six  unknowns:  -P(t),  y{t),  mL{t),  mu{t),  u(0,  t),  and  u{W,  t).  Four  ODE’s  for 
the  first  four  unknowns  are  coupled  with  the  two  functional  equations  (8)  and  (9).  The 
resulting  system  can  be  called  a differential-functional  equation  (DFE)  system.  The 
gradient  matching  method  is  a procedure  for  finding  an  approximate  solution  to  this 
DFE  system  using  a DAE  solver.  It  is  based  on  the  following  semi-group  property  of 
the  heat  equation.  The  temperature  profile  obtained  by  advancing  the  initial  profile 
to  time  t At  IS  identical  with  the  profile  obtained  by  first  advancing  the  initial 
profile  to  time  t and  then  advancing  it  from  time  t to  time  t-\-  At.  Implementation  of 
the  gradient  matching  method  requires  that  storage  be  allocated  for  the  temperature 
profile  at  the  previous  time,  f,  and  at  the  next  time,  t + At.  Given  the  profile  at  time 
t and  values  for  the  six  unknowns  at  time  t + At  (initial  guess  by  the  solver),  the 
profile  G is  advanced  from  time  t to  time  t -f  At.  The  gradient  of  G at  a:  = 0 and 
a:  = VF  is  computed  followed  by  the  residuals  for  the  six  equations  including  equations 
(8)  and  (9).  The  DAE  solver  adjusts  the  stepsize  until  the  residuals  for  all  six 
equations  are  below  an  error  tolerance.  Once  the  solver  has  completed  the  step,  the 
array  storing  the  temperature  profile  for  the  previous  time  is  updated,  and  the  DAE 
solver  is  ready  to  take  its  next  step. 


3.3  Implementation  Details 

The  gradient  matching  method  was  implemented  in  CONRAD2.  For  C0NRAD2  the 
set  of  unknowns  are  constructed  as  follows.  For  each  room  in  which  conduction  is 
modeled,  there  are  four  wall  segments.  The  floor  and  ceiling  are  each  treated  as  a 
segment.  The  vertical  walls  in  the  room  are  divided  into  two  segments,  one  above 
the  layer  interface  and  one  below.  Each  segment  contributes  two  unknowns,  one  for 
the  segment  interior  surface  temperature  and  one  for  the  segment  exterior  surface 
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temperature.  The  four  gas  solution  variables  P,  y,  mu,  and  mi  bring  the  total 
number  of  unknowns  to  twelve  per  room. 

To  compute  the  first  terms  in  equations  (8)  and  (9)  at  time  f + Af,  the  temperature 
profile  must  be  advanced  from  time  t to  time  t At  and  its  endpoint  gradients 
approximated.  To  advance  the  temperature  profile  a simple  MOL  approach  was  used 
in  C0NRAD2.  A graded  (non-uniform)  mesh  with  Uj,  breakpoints  was  introduced  for 
the  spatial  variable  x.  The  second  spatial  derivative  in  the  heat  equation  was  replaced 
by  a second  divided  (finite)  difference  approximation.  This  produces  a system  of  rix— 2 
ode’s  for  the  nx  — 2 unknown  temperatures  at  the  interior  breakpoints.  This  system 
was  solved  by  one  step  of  the  backward  Euler  method.  Crank-Nicholson  was  also 
tried  but  produced  spurious  oscillations  in  the  temperature  profiles  at  the  beginning 
of  the  simulation;  backward  Euler  does  not  suffer  from  this  defect  which  is  related 
to  the  non-uniform  mesh  being  used.  The  solution  at  time  t At  can  be  found  by 
solving  a tridiagonal  system  of  linear  equations.  The  temperature  gradient  at  a:  = 0 
and  time  t + At  was  approximated  by  interpolating  the  temperature  values  at  the 
first  three  breakpoints  with  a quadratic  and  then  evaluating  the  derivative  of  this 
polynomial  at  a:  = 0.  A similar  procedure  was  used  to  approximate  the  temperature 
gradient  at  x = VE  and  time  t At. 

A graded  mesh  scheme  was  chosen  to  allow  breakpoints  to  cluster  near  the  in- 
terior and  exterior  wall  segment  surfaces.  This  is  where  the  temperature  gradi- 
ents are  the  steepest.  A breakpoint  Xf,  was  defined  by  Xb  = min(xp,  kE/2),  where 
Xp  = 2 \/at final  erf c“ ^ ( . 05 ) and  erfc“^  denotes  the  inverse  of  the  complementary  error 
function.  The  value  Xp  is  the  location  in  a semi-infinite  wall  where  the  temperature 
rise  is  5%  after  fftnai  seconds;  it  is  sometimes  called  the  penetration  depth.  Eighty 
percent  of  the  breakpoints  were  placed  on  the  interior  side  of  Xb  and  the  remaining 
twenty  percent  were  placed  on  the  exterior  side. 

The  temperature  profile  in  a conduction  node  at  time  f + Af  is  found  from  the 
temperature  profile  at  time  t by  solving  a tridiagonal  system  as  follows.  Let  u,(f) 
denote  the  approximate  temperature  at  time  t and  breakpoint  Xi,  let  Axj  = x,+i  — Xj 
denote  the  breakpoint  spacing,  and  let  s = The  tridiagonal  system  that  must 

be  solved  is 


1 0 

62  ^2  C2 


^Tix— 1 ^ni  — 1 ('Tii  — 'i 


ui[t  -f  At) 

/o(f  + At) 

U2{t  -f  At) 

U2{t) 

Wnx-1 

^rix  — 1 (0 

Un^(f-t-Af)  ^ 

^ /i(f  + At) 
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where  for  i = 2, . . . , ria;  ~ 1 


* Axi_iAa;,- 


Ax,_i(Ax,_i  + Axi) 

— 5 

^ Ax,(Ax,_i  + Ax,) 

If  the  spacing  is  uniform  (Ax,  = Ax)  and  s = Ax^  then  the  above  coefficient  matrix 
reduces  to 


0 

2 


1 

2 


\ 


V 


1 

2 
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In  wall  materials  composed  of  multiple  slabs,  a different  form  for  a,,  6,  and  c,  must 
be  used  at  the  breakpoints  where  the  slab  material  properties  change.  This  occurs 
because  the  second  partial  of  temperature,  |^,  does  not  exist  at  these  breakpoints. 
Formulas  for  these  coefficients  are  derived  in  Appendix  B using  the  fact  that  heat 
flow  through  a wall  is  continuous.  These  breakpoint  coefficients  are 


a. 


h. 


Ci 

di 


i<r  K- 

Axi  Ax,_i 

K 

Ax,_i 

Ax,- 

0 


where  K^{K~ ) is  the  thermal  conductivity  on  the  right(left)  side  of  breakpoint  x,. 

The  gradients  at  time  f + Af  at  the  points  x = 0 and  x = W are  estimated  from 
quadratics  passing  through  the  first  three  and  last  three  temperature  values.  These 
gradients  are  computed  using  Newton  divided  differences.  Let  u(xo),  u(xi)  and  u(x2) 
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denote  the  temperatures  at  the  first  three  breakpoints.  The  wall  temperature  gradient 
at  the  interior  surface  is  estimated  by 

u[xi,X2]  - u[xi,X2,X3]^Xi 

W 


where 


u[x,y] 

u\x,y,z\ 


u{x)  - u{y) 
x-y 

u[x,y]  - u[y,z] 
X — z 


Similarly  the  wall  temperature  gradient  at  the  exterior  surface  is  given  by 


K 


W 


4 Numerical  Results 

In  this  section  the  MOL  approach  used  in  CONRADl  is  compared  numerically  with 
the  gradient  matching  approach  used  in  CONRAD2.  These  numerical  experiments 
were  conducted  on  a Sun  Sparcstation  2.  CONRADl  and  CONRAD2  were  written 
in  Fortran  77  together  with  two  extensions  that  are  available  on  almost  all  Fortran  77 
compilers.  These  extensions  are  the  “IMPLICIT  NONE”  statement  that  forces  the 
typing  of  all  variables  and  the  “INCLUDE”  statement  that  allows  various  header  files 
to  be  read  at  the  beginning  of  a program  unit.  CONRAD2  has  been  ported  to  and 
tested  on  an  Apple  Macintosh  II,  a Sun  Sparcstation  2,  a Silicon  Graphics  4D35,  and 
an  IBM  Rise  6000  Model  320.  CONRADl  has  been  ported  to  and  tested  on  an  Apple 
Macintosh  II,  a Sun  Sparcstation  2,  and  a Silicon  Graphics  4D35.  The  only  machine 
dependent  parts  of  these  codes  are  the  default  unit  number  for  screen  output,  the 
timing  routines,  and  the  floating  point  constants.  CONRADl  and  CONRAD2  are 
fully  documented  via  comment  statements  which  include  porting  instructions. 

The  user  communicates  with  CONRADl  and  CONRAD2  through  nearly  identical 
input  data  files.  In  this  section,  these  codes  are  compared  on  a four  room  test  case 
for  which  the  CONRADl  input  data  file  has  the  following  structure.  The  units  here 
are  meters  and  seconds. 

'RELATIVE  ERROR  TOLERANCES  FOR  P,  Y,  MU,  ML,  NODETEMP,  NODETEMPGRAD ’ 
l.D-6  l.D-6  l.D-6  l.D-6  l.D-3  l.D-3 

'ABSOLUTE  ERROR  TOLERANCES  FOR  P,  Y,  MU,  ML,  NODETEMP,  NODETEMPGRAD' 
l.D-6  l.D-6  l.D-6  l.D-6  O.DO  O.DO 
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'NUMBER  OF  ROOMS’ 

4 

'FIRE  WATTS, DIM: DATUM  TO  FLOOR, XROOM,YROOM,ZROOM,GXFIRE,GYFIRE,GZFIRE’ 
'FOLLOWED  BY  RADIATION, NODES, C0DES(0=NULL,1=C0NCRETE,2=GYPSUM,3=KA0W00L) ’ 


500000. 

0, 

3.64  3.63 

2.45 

0. 

0. 

0 

.TRUE. 

4 

2 

2 2 2 

0. 

0, 

2.43  18.9 

2.43 

0. 

0. 

0 

.TRUE. 

4 

2 

2 2 2 

0. 

0, 

3.64  3.65 

2.45 

0. 

0. 

0 

.TRUE. 

4 

2 

2 2 2 

0. 

0, 

3.52  3.52 

2.43 

0. 

0. 

0 

.TRUE. 

4 

2 

2 2 2 

'NUMBER  OF  VENTS’ 

6 

'AREA,  DIS:  DATUM  TO  TOP, DATUM  TO  BOTTOM, ROOM  NUMBERS  ON  EACH  SIDE’ 


0.18400E+01 

0.20000E+01 

o 

o 

+ 

u 

o 

o 

o 

o 

o 

o 

1 

2 

0.50000E-01 

0.30000E+00 

0.20000E+00 

1 

-1 

0.60000E-01 

0.20000E+01 

o 

o 

o 

o 

o 

o 

tn 

+ 

o 

o 

2 

3 

0.60000E-01 

0.20000E+01 

o 

o 

o 

o 

o 

o 

w 

+ 

o 

o 

2 

4 

0.50000E-01 

0.30000E+00 

0.20000E+00 

3 

-1 

0.50000E-01 

0.30000E+00 

0.20000E+00 

4 

-1 

'NUMBER  OF  FORCED  VENTS’ 

0 

'AREA,  DIS:  DATUM  TO  TOP, DATUM  TO  VENT, FLOW  RATE, ROOM  NUMBERS  ON  EACH  SIDE’ 
'TFINAL,  TPRINT’ 

600.  60. 

'LOGICAL,  TS,  TF:  DIAGNOSTICS  AT  INTERMEDIATE  STEPS  IN  INTERVAL  [TS,  TF] ’ 

.FALSE.  0.0  0.5 

'NUMBER  OF  BREAKPOINTS  NX’ 

5 


Lines  2 and  4 specify  six  error  tolerances.  For  C0NRAD2  the  last  tolerance 
is  left  out  since  wall  segment  temperature  gradients  are  not  solution  variables  in 
C0NRAD2.  Lines  10,  12,  14,  and  16  specify  that  the  radiation  submodel  should  be 
used  in  each  room,  that  four  wall  segments  (nodes)  should  be  used  in  each  room, 
and  that  all  wall  segment  materials  should  be  gypsum.  Lines  20-25  specify  the  vents 
that  interconnect  the  rooms.  Room  number  —1  indicates  the  outside.  The  last  line 
specifies  the  number  of  breakpoints  used  in  the  wall  segments.  This  number  is  5 
for  CONRAD  1 and  20  for  CONRAD2.  These  numbers  were  chosen  on  the  basis 
of  a set  of  numerical  experiments.  Using  more  than  this  number  of  breakpoints 
does  not  significantly  improve  the  accuracy  of  the  wall  segment  temperature  profile 
computation. 

Table  1 gives  cpu  times  in  seconds  for  three  cases:  gas  equations  only,  gas  plus 
conduction  equations,  gas  plus  conduction  plus  radiation.  The  wall  segment  surface 
temperatures  for  these  two  codes  agreed  to  within  the  error  tolerance  set  in  the  input 
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Table  1:  CPU  Time  Comparisons 


Code 

gas 

gas  -f  cond 

gas  -f-  cond  + rad 

CONRADl 

3.1 

45.2 

97.5 

CONRAD2 

3.4 

30.5 

51.6 

data  file  which  here  is  0.1%.  This  table  illustrates  the  savings  in  cpu  time  provided 
by  CONRAD2  without  any  loss  of  accuracy  in  the  computation  of  wall  segment 
temperatures.  The  number  of  DAE’s  solved  by  DASSL  in  CONRAD  1 is  16  for  the 
gas  case  and  176  for  the  gas  plus  conduction  and  gas  plus  conduction  plus  radiation 
cases.  The  number  of  DAE’s  solved  by  DASSL  in  CONRAD2  is  16  for  the  gas  case 
and  48  for  the  gas  plus  conduction  and  gas  plus  conduction  plus  radiation  cases. 
Notice  that  the  cost  of  using  the  radiation  submodel  is  significant. 

No  attempt  has  been  made  in  CONRAD  1 to  exploit  any  special  structure  in  the 
linear  systems  that  DASSL  solves  at  each  time  step.  DASSL  solves  a nonlinear  system 
at  each  time  step  by  a variant  of  Newton’s  method.  Each  iteration  of  this  method 
requires  the  solution  of  a linear  system.  The  structure  of  the  coefficient  matrix  of 
this  linear  system  is  dependent  on  the  interconnectivity  of  the  rooms.  Consider  an 
example  in  which  four  rooms  are  side  by  side  with  a vent  connecting  rooms  1 and  2,  a 
vent  connecting  rooms  2 and  3,  a vent  connecting  rooms  3 and  4,  and  a vent  from  room 
1 to  the  outside.  In  this  case,  the  coefficient  matrix  will  have  a band  structure  and  the 
band  solver  provided  with  DASSL  can  be  used.  The  minimum  band  width  will  occur 
when  the  rooms  are  in  a “chain”  as  in  this  example.  If  each  room  is  interconnected 
with  all  other  rooms,  this  coefficient  matrix  will  be  full.  Consequently,  there  is  no 
easy  way  to  exploit  the  structure  of  this  coefficient  matrix  without  making  a priori 
assumptions  about  the  interconnectivity  of  the  rooms.  On  the  other  hand,  each  wall 
segment  contributes  a tridiagonal  block  to  this  coefficient  matrix.  Perhaps  it  would 
be  possible  to  create  a special  solver  that  would  be  able  to  exploit  this  feature  of  the 
matrix.  If  CONRAD2  is  later  found  to  have  limitations,  then  creating  a special  solver 
to  use  with  DASSL  in  CONRADl  might  be  worthwhile. 


5 Future  Work 

The  use  of  a DAE  solver  to  solve  a DEE  system  appears  to  be  a new  and  useful  devel- 
opment in  the  application  of  mathematics.  The  theoretical  basis  and  the  limitations 
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of  this  method  will  be  further  explored.  This  is  a mathematical  analysis  issue,  not  a 
programming  issue. 

The  subroutine  CNDUCT  in  C0NRAD2  which  solves  the  heat  equation  for  a wall 
segment  can  be  rewritten  based  on  a piecewise  cubic  Hermite  expansion  approach 
instead  of  the  current  finite  difference  approximation  approach.  This  change  has 
several  advantages.  The  order  of  accuracy  of  the  spatial  approximation  is  increased, 
and  consequently,  fewer  breakpoints  and  less  memory  are  required.  The  expansion 
approach  provides  a natural  method  for  interpolating  the  temperature  profile  between 
breakpoints,  and  it  provides  the  gradients  at  the  endpoints  directly  without  the  need 
for  an  additional  approximation.  The  expansion  approach  should  also  be  less  sensitive 
to  poorly  chosen  breakpoint  spacing.  The  resulting  linear  system  that  must  be  solved 
is  now  pentadiagonal  instead  of  tridiagonal,  but  the  size  of  the  system  is  cut  in  half. 

The  breakpoint  heuristic  in  C0NRAD2  can  be  improved.  The  current  heuristic 
is  based  on  the  final  time  tfinai-  As  an  alternative,  three  different  sets  of  breakpoints 
could  be  generated  based  on  three  times  fshort  < ^mid  < ^fmai-  Subroutine  CNDUCT 
could  be  rewritten  to  use  the  set  based  on  fshort  for  0 < f < fshort,  the  set  based  on  fmid 
for  fshort  < t < tmid,  and  the  set  based  on  tfinai  for  tmid  < t ^ ^finai-  This  would  amount 
to  a crude  moving  mesh  strategy.  When  the  set  of  breakpoints  is  changed,  the  current 
temperatures  at  the  new  breakpoints  must  be  computed.  This  is  especially  easy  with 
the  expansion  approach  since  it  provides  a natural  interpolation  formula. 

Currently,  the  ceiling  and  floor  in  a room  are  each  treated  as  wall  segments  and 
the  vertical  walls  are  divided  into  two  segments,  one  above  the  layer  interface  and  one 
below  the  layer  interface.  Consequently,  the  two  segments  of  the  vertical  walls  have 
surface  areas  which  change  with  time.  The  temperature  distribution  in  the  vertical 
walls  can  be  more  accurately  modeled  by  using  more  than  two  segments  with  fixed 
surface  areas,  or  by  moving  to  a 2-d  conduction  model  in  the  vertical  walls.  Currently, 
the  exterior  surfaces  of  the  wall  segments  in  CONRADl  and  CONRAD2  exchange 
heat  to  ambient.  The  model  could  be  improved  by  allowing  neighboring  rooms  to 
exchange  heat  via  conduction  through  ceilings,  walls,  and  floors.  A simple  first  step 
would  be  to  implement  this  for  the  room(s)  of  fire  origin. 


6 Summary 

Two  standard  methods,  time  splitting  and  the  method  of  lines,  were  presented  for 
coupling  the  partial  differential  equations  that  arise  when  modeling  heat  conduction 
with  the  ordinary  differential  equations  that  arise  in  zone  fire  modeling.  Effective  time 
splitting  implementations  rely  on  the  assumption  that  the  time  split  phenomenon 
changes  over  a significantly  longer  time  scale  than  other  phenomena  of  interest.  This 
assumption  breaks  down  when  wall  materials  are  thin  and/or  highly  conductive  since 
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they  react  more  quickly  to  changing  boundary  conditions.  The  method  of  lines  does 
not  have  these  problems.  However,  efficient  implementations  are  difficult  to  attain 
due  to  the  large  number  of  extra  ordinary  differential  equations  that  are  introduced. 

A third  heat  conduction  coupling  method,  gradient  matching,  was  introduced  to 
address  the  above  deficiencies.  This  method  couples  only  the  surface  wall  tempera- 
tures with  the  other  gas  solution  variables.  These  are  the  temperatures  that  are  of 
primary  concern  since  only  the  wall  surface  temperatures  (not  the  interior  wall  tem- 
peratures) interact  with  the  gas  layers  [via  convective  and  radiative  heat  transfer). 

Numerical  experiments  were  performed  that  demonstrated  that  the  new  scheme 
was  not  only  more  efficient  (as  implemented  in  CONRAD2)  but  agreed  with  estab- 
lished numerical  methods  (using  the  method  of  lines  as  implemented  in  CONRADl) 
to  within  the  numerical  solver  tolerance  error  of  .1  per  cent. 
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Nomenclature 


c 

K 

muimi) 

P 

q" 

t 

Tamb 

i final 

u{x, t) 

X 

Xp 

y 

a 

At 

Axi 

P 


specific  heat  wall  segment  material 

thermal  conductivity  of  wall  segment  material 

mass  of  upper  (lower)  layer 

room  pressure 

net  heat  flux 

time 

ambient  temperature 

final  simulation  time 

temperature  in  a wall  segment 

position  in  wall  segment 

penetration  depth 

layer  interface  height 

thermal  diffusivity,  a = ^ 

time  stepsize 

breakpoint  spacing 

density  of  wall  segment  material 
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A Modeling  Heat  Conduction  Using  the  Method 
of  Lines 


The  standard  MOL  approach  for  solving  the  heat  equation  is  to  discretize  the  spatial 
variable  by  either  replacing  the  spatial  derivatives  with  finite  difference  approxima- 
tions, or  by  expanding  the  unknown  function  as  a linear  combination  of  spatial  basis 
functions  with  time  dependent  coefficients  and  deriving  the  ODE’s  via  collocation. 
CONRAD  1 uses  the  second  approach  and  CONRAD2,  the  first.  The  MOL  generally 
produces  a stiff  system  of  ODE’s.  Although  this  method  was  proposed  many  years 
before,  it  was  not  until  the  advent  of  stiff  ODE  solvers  in  the  1970’s  (see  Gear  [12]) 
that  implementation  of  this  method  was  practical. 

For  each  wall  segment  there  is  an  initial-boundary  value  problem  which  is  coupled 
with  the  gas  ODE’s.  These  equations  are 


du{x,  t) 


dt 


u(x,  0) 


d^u{x, t) 
^ dx^ 


0 < a:  < VF,  f > 0 


Tamb  , ^ < X <W 


9conv.c(0.<)  + 9ild(0.*) 


(10) 

(11) 

(12) 

(13) 


where  u{x^t)  denotes  the  temperature  at  a distance  x into  the  wall  segment  at  time 
t and  Tamb  denotes  the  ambient  temperature.  The  terms  g'onvec  9rld  denote 
the  convective  and  radiative  flux  into  the  interior  and  exterior  surfaces  of  the  wall 
segments.  Finally,  a = — , where  A',  c,  and  p,  denote  the  thermal  conductivity, 
the  specific  heat,  and  the  density  of  the  segment  material.  The  four  gas  ODE’s  are 
coupled  to  the  above  initial-boundary  value  problem  by  the  dependencies  of  ^convec 
and  on  the  gas  layer  temperatures  Tl  and  Tu  which  can  be  determined  from  the 
gas  solution  variables.  In  addition,  the  energy  transfer  rates  qi  and  qu  to  the  lower 
and  upper  gas  layers  contain  terms  to  account  for  the  transfer  of  energy  from  the  gas 
layers  to  the  surfaces  of  the  wall  seqments  via  convection  and  radiation. 


A.l  The  Method  of  Lines  Using  the  Piecewise  Cubic  Her- 
mite  Collocation  Approach,  CONRADl 

It  is  well-known  in  approximation  theory  that  a function  possessing  four  continuous 
derivatives  on  a closed,  finite  interval  can  be  approximated  to  fourth  order  accuracy 
using  piecewise  cubic  Hermite  interpolation.  This  method  of  interpolation  matches 
the  function  and  its  first  derivative  at  a set  of  breakpoints.  Between  the  breakpoints. 
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Figure  3:  Plot  of  Cubic  Hermite  Basis  Functions 


the  function  is  approximated  by  a cubic  polynomial.  The  resulting  approximation 
has  a continuous  first  derivative  (see  [13]).  Let  and  ipi  denote  the  standard  basis 
functions  for  piecewise  cubic  Hermite  interpolation  with  breakpoints  < ...  < Xn^. 
These  basis  functions  have  the  defining  properties 

= 0 

where  Sij  is  the  kroneker  delta  function  with  value  zero  unless  i = j in  which  case  the 
value  is  one. 

For  the  case  where  xi  = — 1,  0:2  = 0,  and  X3  = I,  the  basis  functions  02  and  02 
are  plotted  in  Figure  3 and  are  given  by 

(f)2{x)  = (la:|  - 1)2(21x1  + 1) 

02(a:)  = a;(|x|  - 1)2 

The  MOL  equations  for  CONRADl  are  derived  by  assuming  that  the  wall  segment 
temperature  profiles  have  the  form 

u{x,t)  = ^[ai(f)0,(x)  + 6i(f)0i(x)]  . (14) 

i=l 
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Equation  (14)  is  then  substituted  into  equation  (10)  to  obtain 


Tlx 

S (a'(t)(^.(x)  + h'i{t)xl)i{x))  = ocJ^  {ai{t)(t>';{x)  + hi{t)xl)';{x))  . (15) 

i i 

Next,  2nx  — 2 ODE’s  involving  the  2nx  unknown  coefficients,  a,(f),  6,(f),  f = 1, . . . , Ua,, 
are  derived  by  requiring  that  equation  (15)  be  satisfied  at  the  following  two  Gaussian 
points 


P2J-1  = Xj  + 5— - xj), 
3- v^, 

P2j  ^j+1  ^ \^j+i  ^j) 


in  each  subinterval  j = 1,  • • • , rix  — 1.  The  boundary  conditions  (12)  and  (13)  provide 
two  additional  algebraic  equations  given  by 


-/C[a,(()«i',(0)  + 6i(()V’!(0)]  = 9conv,c(0.*)  + 9"ad(0.<). 

These  discretizations  lead  to  a system  of  2ux  — 2 ODE’s  and  two  algebraic  equations 
of  the  form 


= aBy  + F, 


where  A and  B are  2nx  x 2nx  matrices  and  y and  F are  2nx  vectors.  Analytic 
differentiation  of  the  boundary  conditions  with  respect  to  time,  which  would  lead 
to  a system  of  2nx  ODE’s,  is  not  practical  because  the  radiation  flux  terms  in  the 
boundary  conditions  are  implicitly  defined,  nonlinear  function  of  the  gas  variables,  P, 
y,  mui  and  mi,  and  the  wall  segment  surface  temperatures.  As  a result  the  boundary 
conditions  must  be  treated  as  algebraic  equations.  For  = 4 and  q”  = q"onvec  + 
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the  vectors  y and  F have  the  form 


/ \ 
ai 

02 


03 

^3 

0,4 

< / 


and  the  matrices  A and  B have  the  form 


^ 0 

0 

0 

0 

0 

0 

0 

0 ^ 

</»2(Pl) 

V’2(Pl) 

0 

0 

0 

0 

<t>\{V2) 

V’l(P2) 

(t>2{P2) 

V’2(P2) 

0 

0 

0 

0 

0 

0 

^2(^3) 

V’2(P3) 

Mps) 

Mp3) 

0 

0 

0 

0 

4>2[P4) 

^2{Pa) 

Mp4) 

Mp4) 

0 

0 

0 

0 

0 

0 

Mps) 

Mps) 

Mps) 

Mp3) 

0 

0 

0 

0 

Mpe) 

Mpe) 

Mpe) 

Mpe) 

V 0 

0 

0 

0 

0 

0 

0 

0 y 

{ q"{0,t)  \ 

pc 

0 
0 
0 
0 
0 
0 

pc 


(16) 
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and 


' 0 

1 

0 

0 

0 

0 

0 

0 ^ 

<f>'^{pl) 

0 

0 

0 

0 

^"{P2) 

<I>2{P2) 

V’2(P2) 

0 

0 

0 

0 

0 

0 

4>'i{p.) 

^2{P3) 

<f>3{P3) 

^3{P3) 

0 

0 

0 

0 

<t>2ip^) 

^2{P4) 

<^3(P4) 

V^3(P4) 

0 

0 

0 

0 

0 

0 

<?^3(P5) 

V^sIPs) 

<^4(P5) 

V’^'Ips) 

0 

0 

0 

0 

(f>3{P6) 

V’sCPe) 

^':{pe) 

^(pe) 

V 0 

0 

0 

0 

0 

0 

0 

1 > 

The  DAE  solver  DASSL  by  Petzold  [7,  8]  was  chosen  for  use  in  CONRAD  1 and 
CONRAD2.  It  is  the  most  widely  used  production  code  for  DAE’s  at  this  time.  The 
following  brief  description  of  DASSL  is  taken  from  [7].  DASSL  is  a code  for  solving 
systems  of  DAE’s  of  the  form 


= 0 

(18) 

y{to) 

= yo 

(19) 

y'{to) 

= y'o^ 

(20) 

where  F,  y,  and  y'  are  A^-dimensional  vectors.  The  basic  idea  for  solving  DAE  systems 
using  numerical  ODE  methods  is  to  replace  the  derivative  in  (18)  by  a difference 
approximation,  and  then  to  solve  the  resulting  system  for  the  solution  at  the  current 
time  fn+i  using  Newton’s  method.  For  example,  replacing  the  derivative  in  (18)  by 
the  first  order  backward  difference,  we  obtain  the  implicit  Euler  formula 

f ~^")  =0,  (21) 

V «n+l  / 

where  /i^+i  = ^n+i  ~ in-  This  nonlinear  system  is  then  usually  solved  using  some 
variant  of  Newton’s  method.  The  algorithms  used  in  DASSL  are  an  extension  of  this 
basic  idea.  Instead  of  always  using  the  first  order  formula  (21),  DASSL  approximates 
the  derivative  using  the  A:th  order  backward  differentiation  formula,  where  k ranges 
from  one  to  five.  At  every  step  it  chooses  the  order  k and  the  stepsize  hn+i,  based 
on  the  behavior  of  the  solution.  DASSL  can  solve  index  zero  and  one  systems.  The 
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Figure  4:  Ceiling  Temperature  Profiles  at  20,  60  and  120  Seconds 


index  of  the  DAE  system  (18)  is  the  minimum  number  of  times  that  all  or  part  of 
this  system  must  be  differentiated  with  respect  to  t in  order  to  determine  y'  as  a 
continuous  function  of  y and  t. 


A. 2 Graded  Spatial  Meshes 

The  MOL  chooses  the  time  discretization  to  maintain  accuracy  and  stability,  but 
the  user  must  choose  the  spatial  discretization.  During  the  first  seconds  of  a fire 
simulation,  wall  segment  temperature  profiles  typically  have  steep  gradients  near  the 
interior  wall  segment  surfaces.  Consequently,  uniform  meshes  are  not  efficient.  For 
both  CONRAD  1 and  CONRAD2,  graded  meshes  were  developed  with  the  grading 
dependent  on  the  final  simulation  time  ffinai-  These  graded  mesh  were  optimized  for 
the  case  when  the  fire  energy  release  rate  takes  a step  jump  t = 0 and  then  is 
constant  thereafter.  In  this  case,  the  steepest  temperature  gradients  occur  near  the 
interior  surfaces  of  the  wall  segments.  As  the  simulation  evolves,  these  temperature 
profiles  tend  to  flatten  out.  Figure  4 shows  the  ceiling  temperature  profiles  at  various 
times  during  a one  room  simulation.  The  room  is  3 m long,  2 m wide,  3 m high 
and  contains  a 1 Mw  fire  on  the  floor.  It  has  a single  1 m?  vent  to  the  outside.  All 
conduction  nodes  are  made  from  gypsum.  The  four- wall  radiation  model  documented 
in  [14]  was  used  to  model  radiation  heat  transfer.  The  simulation  runs  for  2 minutes. 
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The  general  qualitative  features  of  the  profiles  in  Figure  4 are  exhibited  by  the 
semi-infinite  (0  < a:  < oo)  rod  solution  to  the  heat  equation 


— ^amb  d”  (^1  Tajjjt))Grfc 


(22) 


where  u{0,t)  = uj,  a constant,  the  initial  temperature  is  T^mh,  ^ erfc(x)  = 

e~^^dt  denotes  the  complementary  error  function.  Using  this  solution  as  a 
guide,  a graded  mesh  was  developed  to  accommodate  both  short  times  and  times 
near  f final  using  Ux  > 3 breakpoints  for  CONRAD  1 and  > 10  breakpoints  for 
CONRAD2.  In  both  cases  a:i  = 0 and  a:^  = It  is  impossible  for  a fixed  mesh 
to  be  optimum  for  all  times  0 < f < ffinai-  The  meshes  used  in  CONRAD  1 and 
CONRAD2  were  designed  to  concentrate  most  of  the  breakpoints  between  0 and  a 
breakpoint 

W 

Xb  = mm(a:p,  y)  . 

The  point  Xb  was  chosen  to  be  the  smaller  of  the  midpoint  of  the  wall  segment,  VF/2, 
and  the  penetration  depth,  Xp,  defined  by 

Xp  :=  2\/Q!^fmai erfc“^(.05)  . 


Table  2 gives  penetration  depths  for  several  wall  materials  with  ffinai  = 600s.  To 
obtain  penetration  depths  for  a different  ffinai,  multiply  the  values  in  column  five  of 
Table  2 by  final- 

In  CONRAD2,  80  percent  of  the  breakpoints  are  to  the  left  of  Xb-  The  breakpoints 
on  either  side  of  Xb  are  quadratically  graded  so  that  they  cluster  near  x = 0 and 
X = VF. 

The  breakpoint  design  for  CONRAD2  is  somewhat  different.  For  a fixed  time  t > 
0,  the  heat  equation  solution  in  (22)  is  flattened  out  for  x > X{t)  = 2\/^erfc“^(.05). 
To  resolve  the  temperature  profile  at  this  time,  a breakpoint  should  be  placed  near 
X{t).  Two  of  the  breakpoints  are  generated  this  way.  The  breakpoint  X2  is  the 
minimum  of  X(fprint)  and  ^ with  the  ^ term  required  to  provide  short  time  accuracy 
when  fprint  is  not  sufficiently  small.  CONRAD  I and  CONRAD2  print  out  data  every 
Sprint  seconds.  The  breakpoint  Xn,_i  is  set  to  xj,.  By  the  time  X{t)  has  reached 
the  ceiling  midpoint,  the  temperature  profile  for  x > ^ Is  linear  enough  so  that 
breakpoints  at  ^ ^^d  W suffice. 
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Table  2:  Penetration  Depths,  tfinai  = 6OO5 


Material 

K 

{WI{mK)) 

P 

[kglm^) 

c 

{WsHkgK)) 

2y/aeric  ^(.05) 

(mV^) 

Xp 

Copper 

387 

380 

8940 

2.96  X 10-2 

0.72 

Oak 

.17 

2380 

800 

8.28  X lO-'* 

2.03  X 10-2 

PMMA 

.19 

1420 

1190 

9.29  X 10“^ 

2.28  X 10-2 

Brick 

.69 

840 

1600 

1.99  X 10-3 

4.86  X 10-2 

Kaowool 

.22 

128 

1047 

3.55  X 10-3 

8.70  X 10-2 

Gypsum 

.16 

800 

900 

1.31  X 10-3 

3.20  X 10-2 

Concrete 

1.75 

2200 

1000 

2.47  X 10-3 

6.06  X 10-2 

B Heat  Conduction  in  Multiple  Slab  Wall  Seg- 
ments 

In  this  section  methods  are  presented  for  handling  the  case  of  a wall  segment  composed 
of  slabs  of  different  materials.  For  simplicity,  a two  slab  wall  segment  is  examined. 

B.l  The  Heat  Equation 

For  a two  slab  wall  segment,  the  initial-boundary  value  problem  defined  in  equations 
(3)-(6)  must  be  replaced  by  the  following  interface  problem.  Find  u(x,f)  so  that 


du{x,  t) 
dt 

= 0<x<L<  VF,  f>0 

ox^ 

(23) 

du{x^  t) 
dt 

OX^ 

(24) 

u{x^  0) 

= Tamb,  0 < X < IF 

(25) 

u(0, t) 

= kit),  t > 0 

(26) 

u[W,t) 

= Mt),  t > 0 

(27) 

u{L—,t) 

(28) 

du~ 

dx 

(29) 
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Here,  u{L—,t)  and  u{L-\-^t)  denote  left  and  right  hand  limits,  while  and 
denote  left  and  right  hand  derivatives.  Material  1 occupies  the  interval  [0,L],  while 
material  2 occupies  the  interval  [L,W].  The  interface  conditions  (28)  and  (29)  re- 
quired continuity  of  temperature  and  heat  flux  across  the  interface  at  x = L.  Equation 
(29)  assumes  that  there  is  no  heat  source  embedded  in  the  interface;  this  equation  is 
derived  from  conservation  of  energy. 


B.2  Finite  Difference  Equations 


Introduce  breakpoints  0 = Xj  < . . . < x^  = W,  with  = L for  some  m,  1 < m < n, 
and  let  Axj  = Xj^i  —Xj.  Let  the  temperature  at  breakpoint  Xj  at  time  t be  denoted  by 
Uj{t).  Replacing  the  second  spatial  derivative  in  the  heat  equations  (23)  and  (24)  by 
a second  divided  difference  approximation  and  the  first  spatial  derivative  in  equation 
(29)  by  a first  divided  difference  approximation,  yields  the  following  DAE  system: 


/ 


Uj(0) 

Uo 


Un 


^m— 1 


a 


(1) 


a 


(2) 


“j+l— “J 

7 

1 

1 

Axj_i 

Ax,+Ax,_i 

2 

Uj-Uj-i 

Axj 

Axj_i 

Ax,-|-Ax,_i 

Tamb  1 J — 1 , . . . , n 
/o 


i = 2, . . . , m - 1 
j = m-|-l,...,n— 1 


/i 

_/^(2) 


^m+1 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


If  the  backward  Euler  method  (see  equation  (21))  is  used  to  advance  the  solution 
from  time  t to  time  t 4-  Af,  the  following  tridiagonal  linear  system  arises: 


bjUj—\(t  "b  Ai)  T (ijUj{t  At)  CjUj^i(t  T At)  — dj  . (^h) 


The  coefficients  are  as  follows.  For  j = 2, . . . , m — 1, 


Qj 


dj 


Axj_i(Axj  T Axj_i) 
AxjAxj_i 


Axj(Axj  + Axj-i) 
Uj{t)  . 
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For  j = m+  l,...,n  — 1, 


Axj_i(Axj  + Axj_i) 

a(2) 

1 + ~7 T 

-2a(") 

Axj{Axj  + Axj^i) 

Uj{t)  ■ 

For  j = m, 

Axj_i 

Axj_i  Axj 

Axj 

0 . 

B.3  Cubic  Hermite  Expansion  Equations 

An  expansion  of  the  form  (14)  can  be  used  for  each  slab.  Introduce  breakpoints 
0 = Xi^\...,x^^  = L and  coefficients  for  the  interval  [0,  L],  and  introduce 

breakpoints  L = . , . , = W and  coefficients  for  the  interval  [L,  W]. 

Two  DAE  systems  arise  from  applying  the  cubic  Hermite  collocation  procedure  of 
Section  A.l  to  equations  (23),  (25),  and  (26),  and  to  equations  (24),  (25),  and  (27). 
These  DAE  systems  are  coupled  by  the  interface  conditions  (28)  and  (29)  which 
discretize  as 
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